The present invention relates to radioactive-emission measurements of a volume. More particularly, the present invention relates to the accurate reconstruction of the volume, based on measurements from non-uniform views of a volume, and on the dynamic selection of views during the data acquisition process. Of particular interest is view selection for medical imaging and/or in conjunction with medical instruments, such as guided minimally invasive surgical instruments.
Radionuclide imaging is one of the most important applications of radioactivity in medicine. Its purpose is to obtain a distribution image of a radioactively labeled substance, e.g., a radiopharmaceutical, within the body following administration thereof to a patient. Radioactive-emission imaging relies on the fact that in general, pathologies, such as malignant tumors, malfunctioning organs, and inflammations, display a level of activity different from that of healthy tissue. Thus, radiopharmaceuticals, which circulate in the blood stream, are picked up by the active pathologies to a different extent than by the surrounding healthy tissue; in consequence, the pathologies are operative as radioactive-emission sources and may be detected by radioactive-emission imaging. It will be appreciated that the pathology may appear as a concentrated source of high radiation, a hot region, as may be associated with a tumor, or as a region of low-level radiation, which is nonetheless above the background level, as may be associated with carcinoma.
A reversed situation is similarly possible. Dead tissue has practically no pick up of radiopharmaceuticals, and is thus operative as a cold region.
Thus radiopharmaceuticals may be used for identifying active pathologies as well as dead tissue.
In the discussion that follows, the term body structure is intended to include organs, portions of organs, a part of the body, and a whole body. The term organ target is intended to include pathological features within organs. These pathological features may be expressed, by radioactive-emission imaging, as any one of the following:
i. hot regions, of a radioactive emission intensity higher than the background level;
ii. regions of low-level radioactive emission intensity, which is nonetheless above the background level; and
iii cold regions, of a radioactive emission intensity, lower than the background level.
Examples of radiopharmaceuticals include monoclonal antibodies or other agents, e.g., fibrinogen or fluorodeoxyglucose, tagged with a radioactive isotope, e.g., 99Mtechnetium, 67gallium, 201thallium, 111indium, 123iodine, 125iodine and 18fluorine, which may be administered orally or intravenously. The radiopharmaceuticals are designed to concentrate in the area of a tumor, and the uptake of such radiopharmaceuticals in the active part of a tumor, or other pathologies such as an inflammation, is higher and more rapid than in the tissue that neighbors the tumor. Thereafter, a radiation-emission-measuring-probe, which may be configured for extracorporeal or intracorporeal use, is employed for locating the position of the active area. Another application is the detection of blood clots with radiopharmaceuticals such as ACUTECT from Nycomed Amersham for the detection of newly formed thrombosis in veins, or clots in arteries of the heart or brain, in an emergency or operating room. Yet other applications include radioimaging of myocardial infarct using agents such as radioactive anti-myosin antibodies, radioimaging specific cell types using radioactively tagged molecules (also known as molecular imaging), etc.
The usual preferred emission for such applications is that of gamma rays, which emission is in the energy range of approximately 11-511 KeV. Beta radiation and positrons may also be detected.
Radioactive-emission imaging is performed with a radioactive-emission-measuring detector, such as a room temperature, solid-state CdZnTe (CZT) detector, which is among the more promising that are currently available. It may be configured as a single-pixel or a multi-pixel detector. Alternatively, another solid-state detector such as CdTe, HgI, Si, Ge, or the like, or a combination of a scintillation detector (such as NaT(Tl), LSO, GSO, CsI, CaF, or the like) and a photomultiplier, or another detector as known, may be used.
1a-1i schematically illustrate detecting units 102 and detecting blocks 101 of various geometries and constructions, and radioactive-emission-measuring probes associated with them.
FIG. 1a schematically illustrates a detecting unit 102, formed as a single-pixel detector 104, for example, a room-temperature solid-state CdZnTe (CZT) detector, having a diameter D and a thickness τd. Both the detector diameter D, or a diameter equivalent in the case of a non-circular detector, and the detector thickness τd affect the detecting efficiency. The detector diameter determines the surface area on which radioactive emission impinges; the greater the surface area, the greater the efficiency. The detector thickness affects the stopping power of the detector. High energy gamma rays may go through a thin detector, and the probability of their detection increases with detector thickness. By itself, a single-pixel detector cannot generate an image; rather, all counts are distributed over the surface area of the detector.
FIG. 1b schematically illustrates the detecting unit 102 with a collimator 108, formed as a single cell of a diameter D, a length L, and a septa thickness τ, attached to the detector 104. The collimator 108 may be, for example, of lead, tungsten or another material which substantially blocks gamma and beta rays.
The collimator's geometry, and specifically, the ratio of D/L, provides the detecting unit 102 with a collection angle δ analogous to a viewing angle of an optical camera. The collection angle δ limits the radioactive-emission detection to substantially only that radioactive emission, which impinges on the detector 104 after passing through a “corridor” of the collimator 108 (although in practice, some high-energy gamma rays may penetrate the collimator's walls).
FIG. 1c schematically illustrates a block 101 of the detecting units 102, with the collimator 108, formed as a multi-cell collimator, of a cell diameter D. The collection angle δ is defined for each of the detecting units 102 in the block, and each of the detecting units 102 forms a pixel in the block 101.
FIG. 1d schematically illustrates a radioactive-emission-measuring probe 100 which comprises several detecting units 102, of different geometries and different collection angles δ, within a housing 107.
FIGS. 1e-1i schematically illustrate the block 101, formed as a combination of a scintillation detector (such as NaI(Tl), LSO, GSO, CsI, CaF, or the like), a collimator grid, and photomultipliers.
As seen in FIG. 1e, the block 101, having proximal and distal ends 109 and 111, respectively, vis a vis an operator (not shown), is formed of the scintillation detector 104, of a single pixel, and the collimators 108, to create the detecting units 102. A plurality of photomultipliers 103 is associated with the single pixel scintillation detector 104, and with proper algorithms, as known, their output can provide a two dimensional image of the scintillations in the single pixel scintillation detector 104. In essence, this is an Anger camera, as known.
The distal view 111 of the collimator grid is seen in FIG. 1f. 
Two optional proximal views 109 of the photomultipliers 103 are seen in FIGS. 1g and 1h, as a square grid arrangement, and as an arrangement of tubes.
An Anger camera 117, of the block 101 in the housing 107 is seen in FIG. 1i. 
In each of the cases of FIGS. 1a-1i, the geometry of the collimator 108 determines the collection angle δ, wherein with no collimator, the collection angle δ, is essentially a solid angle of 4π steradians. Thus, the collimator's geometry affects both the detection efficiency and the image resolution, which are defined as follows:
i. The detection efficiency is the ratio of measured radiation to emitted radiation; and
ii. The image resolution is the capability of making distinguishable closely adjacent radioactive-emission organ targets, or the capability to accurately determine the size and shape of individual radioactive-emission organ targets.
Naturally, it is desired to optimize both the detection efficiency and the image resolution. Yet, they are inversely related to each other. The detection efficiency increases with increasing collimator's collection angle, and the image resolution decreases with increasing collimator's collection angle. For example, when the ratio of D/L is 1/2, the collection angle δ is substantially 2.5 steradians, so the cell views incident radiation within the confinement of about a 2.5-steradian sector. However, when the ratio of D/L is 1/12, the collection angle δ is substantially 0.31 steradians, so the cell views incident radiation within the confinement of about a 0.31-steradian sector.
Once the emission data is obtained, the data is processed to reconstruct the intensity distribution within the measured volume. The reconstruction process is generally complex, due to the large quantity of data which must be processed in order to obtain an accurate reconstruction. The following statistical model may be used to perform reconstruction.
We assume an intensity distribution, I, defined over an input space U, where U comprises a set of basic elements (e.g., pixels in two dimensional spaces, voxels in three dimensional spaces), and I(u) is the intensity of a given basic element uεU. A detecting unit positioned on a radiation-emission-measuring-probe takes a series of measurements y=(yi)i=1T from different positions and orientations around the volume U. The geometrical and physical properties of the detecting unit, together with its position and orientation in a given measurement i, determine the detection probability φi(u) of a photon emitted from location u. Thus the effective intensity of location u as viewed by the detecting unit during measurement i is φi(u)I(u).
The random count Xi(u) of photons that are emitted from location u and detected in measurement i is modeled by a Poisson process with mean φi(i)I(u). The total count of photons detected in measurement i is thus:yi˜Poisson(ΣuεUφi(u)I(u))  (1a)or in matrix notation:y=Poisson(ΦI)  (1b)where y is the vector of measurements yi, Φ is a matrix of detection probabilities over measurements i and voxels u, and I is a vector of intensity per voxel u. The reconstruction problem is to reconstruct the intensities I from the measurements y.
The 2-D Radon transform is a mathematical relationship which may be used to reconstruct the emission intensities of volume U when the set of measurements (yt)t=1T is unconstrained. The Radon transform is not statistical and does not take into account the Poissonian nature of the counts. In addition, it models the views as line projections. The Radon transform maps the spatial domain (x,y) to the Radon domain (p,φ). For a fixed projection angle, the Radon transform is simply a projection of the object. A technique known in the art as filtered back-projection (FBP) uses a back-projection operator and the inverse of the Radon transform to reconstruct the intensity distribution in volume U from measurements (yt)t=1T.
The basic, idealized problem solved by the FBP approach is to reconstruct an image from its Radon transform. The Radon transform, when properly defined, has a well-defined inverse. However, in order to invert the transform one needs measured data spanning 180°. In many medical imaging situations, the positioning of the detecting unit relative to the emitting object is constrained, so that complete measured data is not available. Reconstruction based on filtered back-projection is therefore of limited use for medical imaging. Maximum likelihood (ML) and Maximum A Posteriori (MAP) estimation methods, which address the statistical nature of the counts, have been found to provide better image reconstructions than FBP.
Limited-angle tomography is a reconstruction technique in the related art which reconstructs an image from projections acquired over a limited range of angular directions. The success of the reconstruction process depends upon the extent of the angular range acquired compared with the angular range of the missing projections. Any reconstruction from a limited range of projections potentially results in spatial distortions (artifacts) in the image. Limited angle techniques can be applied for both the Radon transform and the statistical models, but better results are generally achieved within the statistical framework. While it is known that the severity of the artifacts increases with the increasing angular range of the missing projections, limited-angle tomography does not provide information on which projections should be used in order to most effectively reconstruct the image.
Maximum likelihood (ML) estimation is a widely used method in the related art for reconstructing an image from a constrained set of measurements. A parameterization of the generative model described above is obtained by assigning an intensity I(u) to every voxel in U. The likelihood of the observed data y=(yt)t, given the set of parameters I={I(u):uεU} is:
                                                                        L                ⁡                                  (                                      y                    ❘                    I                                    )                                            =                            ⁢                              ln                ⁢                                                                  ⁢                                  P                  ⁡                                      (                                          y                      ❘                      I                                        )                                                                                                                          =                            ⁢                              ln                ⁢                                                      ∏                    t                                    ⁢                                      P                    ⁡                                          (                                                                        y                          t                                                ❘                        I                                            )                                                                                                                                              =                            ⁢                                                ∑                  t                                ⁢                                  ln                  ⁢                                                                          ⁢                                      P                    (                                                                                            ∑                          u                                                ⁢                                                                              x                            t                                                    ⁡                                                      (                            u                            )                                                                                              ❘                      I                                        )                                                                                                                          =                            ⁢                                                ∑                  t                                ⁢                                  ln                  ⁢                                                                          ⁢                                      Poisson                    (                                                                  y                        t                                            ❘                                                                        ∑                          u                                                ⁢                                                                                                            ϕ                              t                                                        ⁡                                                          (                              u                              )                                                                                ⁢                                                      I                            ⁡                                                          (                              u                              )                                                                                                                                            )                                                                                                                          =                            ⁢                                                ∑                  t                                ⁢                                  {                                                            -                                                                        ∑                          u                                                ⁢                                                                                                            ϕ                              y                                                        ⁡                                                          (                              u                              )                                                                                ⁢                                                      I                            ⁡                                                          (                              u                              )                                                                                                                                            +                                                                  y                        t                                            ⁢                      ln                      ⁢                                                                        ∑                          u                                                ⁢                                                                                                            ϕ                              t                                                        ⁡                                                          (                              u                              )                                                                                ⁢                                                      I                            ⁡                                                          (                              u                              )                                                                                                                                            -                                          ln                      ⁡                                              (                                                                              y                            t                                                    ⁢                          1                                                )                                                                              }                                                                                        (        2        )            Note that the lower and upper bound of an indexing variable (such as voxels u and time index t) are omitted in the following description, when they are clear from the context.
There is currently no analytic way to solve Eqn. 2 for the maximum of the likelihood function. However, optimization methods that find local maxima of the likelihood are known. One such method is the Expectation-Maximization (EM) process. In EM estimation, there is no guarantee that the sequence converges to a maximum likelihood estimator. For multimodal distributions, this means that an EM algorithm will converge to a local maximum (or saddle point) of the observed data likelihood function, depending on starting values.
Since the data generated by the model is only partially observable by our measurements, a basic ingredient of the Expectation-Maximization formalism is to define a set of random variables that completely define the data generated by the model. In the current case, since Yt=ΣuXt(u), the set of variables {Xu(t):uεU; t=1, . . . , T} is such a set; the generated data is x=(xt)t where xt=(xt(u))u, and the observed data y is completely determined by x. The main tool in the EM formalism is the complete data likelihood:
                                                                        ln                ⁢                                                                  ⁢                                  P                  ⁡                                      (                                          x                      ❘                      I                                        )                                                              =                            ⁢                              ln                ⁢                                                      ∏                    t                                    ⁢                                      P                    ⁡                                          (                                                                        x                          t                                                ❘                        I                                            )                                                                                                                                              =                            ⁢                                                ∑                  t                                ⁢                                  ln                  ⁢                                                                          ⁢                                                            ∏                      u                                        ⁢                                          Poisson                      ⁡                                              (                                                                                                            x                              t                                                        ⁡                                                          (                              u                              )                                                                                ❘                                                                                                                    ϕ                                t                                                            ⁡                                                              (                                u                                )                                                                                      ⁢                                                          I                              ⁡                                                              (                                u                                )                                                                                                                                    )                                                                                                                                                                    =                            ⁢                                                ∑                  t                                ⁢                                                      ∑                    u                                    ⁢                                      {                                                                                            -                                                                                    ϕ                              t                                                        ⁡                                                          (                              u                              )                                                                                                      ⁢                                                  I                          ⁡                                                      (                            u                            )                                                                                              +                                                                                                    x                            t                                                    ⁡                                                      (                            u                            )                                                                          ⁢                                                  ln                          ⁡                                                      (                                                                                                                            ϕ                                  t                                                                ⁡                                                                  (                                  u                                  )                                                                                            ⁢                                                              I                                ⁡                                                                  (                                  u                                  )                                                                                                                      )                                                                                              +                                              ln                        ⁡                                                  (                                                                                                                    x                                t                                                            ⁡                                                              (                                u                                )                                                                                      ⁢                            1                                                    )                                                                                      }                                                                                                          (        3        )            
Since the likelihood depends on the complete data, which is only partially observable, we take its expectation with respect to the space of the unobserved data, given the current set of hypothesized parameters (i.e. the current estimator). The result is a function Q(I|I′) which assigns likelihood to sets I of model parameters, given the current set I′, and given the observed data y:
                                                                        Q                ⁡                                  (                                      I                    ❘                                          I                      ′                                                        )                                            =                            ⁢                              E                ⁡                                  [                                                                                    ln                        ⁢                                                                                                  ⁢                                                  P                          ⁡                                                      (                                                          x                              ❘                              I                                                        )                                                                                              ❘                      y                                        ;                                          I                      ′                                                        ]                                                                                                        =                            ⁢                                                ∑                  t                                ⁢                                                      ∑                    u                                    ⁢                                      {                                                                                            -                                                                                    ϕ                              t                                                        ⁡                                                          (                              u                              )                                                                                                      ⁢                                                  I                          ⁡                                                      (                            u                            )                                                                                              +                                                                        E                          ⁡                                                      [                                                                                                                                                                x                                    t                                                                    ⁡                                                                      (                                    u                                    )                                                                                                  ❘                                                                  y                                  t                                                                                            ;                                                              I                                ′                                                                                      ]                                                                          ⁢                                                  ln                          ⁡                                                      (                                                                                                                            ϕ                                  t                                                                ⁡                                                                  (                                  u                                  )                                                                                            ⁢                                                              I                                ⁡                                                                  (                                  u                                  )                                                                                                                      )                                                                                              +                      C                                        }                                                                                                          (        4        )            where C is a term which is independent of the intensities I. The function Q(I|I′) is maximized by the following new estimates:
                                          I            ⁡                          (              u              )                                =                                    1                                                ∑                  t                                ⁢                                                      ϕ                    t                                    ⁡                                      (                    u                    )                                                                        ⁢                                          ∑                t                            ⁢                              E                ⁡                                  [                                                                                                              x                          t                                                ⁡                                                  (                          u                          )                                                                    ❘                                              y                        t                                                              ;                                          I                      ′                                                        ]                                                                    ;                  ∀                      u            ∈                          U              .                                                          (        5        )            
The expectation in Eqn. 5 is obtained as follows:
                              ⁢                                    P                                                X                  t                                ⁡                                  (                  u                  )                                                      ⁡                          (                                                                    ⁢                                          (                      u                      )                                                        ❘                                ;                                  I                  ′                                            )                                      =                                            ⁢                              (                                                      ❘                                          ⁢                                              (                        u                        )                                                                              ;                                      I                    ′                                                  )                            ⁢              ⁢                              (                                                      ⁢                                          (                      u                      )                                                        ❘                                      I                    ′                                                  )                                                    ⁡                              (                                  ❘                                      I                    ′                                                  )                                              =                                                                                                                Poisson                      (                                                                                                    y                            t                                                    -                                                                                    x                              t                                                        ⁡                                                          (                              u                              )                                                                                                      ❘                                                                              ∑                                                          v                              ≠                              u                                                                                ⁢                                                                                                          ⁢                                                                                                                    ϕ                                t                                                            ⁡                                                              (                                v                                )                                                                                      ⁢                                                                                          I                                ′                                                            ⁡                                                              (                                v                                )                                                                                                                                                        )                                                                                                                                  Poisson                      ⁡                                              (                                                                              ⁢                                                          (                              u                              )                                                                                ❘                                                                                    ϕ                              t                                                        ⁢                            ⁢                                                          (                              u                              )                                                        ⁢                                                                                          I                                ′                                                            ⁡                                                              (                                u                                )                                                                                                                                    )                                                                                                                        Poisson                ⁡                                  (                                                            y                      t                                        ❘                                          ⁢                                                                        ϕ                          t                                                ⁢                        ⁢                                                  (                          v                          )                                                ⁢                                                                              I                            ′                                                    ⁡                                                      (                            v                            )                                                                                                                                )                                                      =                          Binomial              ⁡                              (                                                                            ⁢                                              (                        u                        )                                                              ❘                                                                  ⁢                                                  (                          u                          )                                                ⁢                                                                              I                            ′                                                    ⁡                                                      (                            u                            )                                                                                                                      ⁢                                                  ϕ                          t                                                ⁢                        ⁢                                                  (                          v                          )                                                ⁢                                                                              I                            ′                                                    ⁡                                                      (                            v                            )                                                                                                                                ;                                      y                    t                                                  )                                                                        (        6        )            
It follows that E[xt(u)|yt;
                    I        ′            ]        =                  y        t            ⁢                                                  ϕ              t                        ⁡                          (              u              )                                ⁢                                    I              ′                        ⁡                          (              u              )                                                            ∑            v                    ⁢                                          ⁢                                                    ϕ                t                            ⁡                              (                v                )                                      ⁢                                          I                ′                            ⁡                              (                v                )                                                          ,and hence the EM iteration is:
                              I          ⁡                      (            u            )                          =                              1                                          ∑                t                            ⁢                                                          ⁢                                                ϕ                  t                                ⁡                                  (                  u                  )                                                              ⁢                                    ∑              t                        ⁢                                          y                t                            ⁢                                                                                          ϕ                      t                                        ⁡                                          (                      u                      )                                                        ⁢                                                            I                      ′                                        ⁡                                          (                      u                      )                                                                                                            ∑                    v                                    ⁢                                                                          ⁢                                                                                    ϕ                        t                                            ⁡                                              (                        v                        )                                                              ⁢                                                                  I                        ′                                            ⁡                                              (                        v                        )                                                                                                                                                    (        7        )            
It is provable that each EM iteration improves the likelihood. Thus, given a random starting estimator, the EM algorithm iterates the above improvement step until it converges to a local maximum of the likelihood. Several random starts increase the chance of finding a globally good estimator.
It is usually desired to maximize the expected posterior probability (given a proper prior) rather than the expected likelihood. In that case we assume a prior probability on the intensities P(I)=ΠuP(I(u)). A proper conjugate prior for the Poisson distribution is the Gamma distribution:
                              P          ⁡                      (                          I              ⁡                              (                u                )                                      )                          =                              Gamma            ⁡                          (                                                I                  ⁡                                      (                    u                    )                                                  ❘                            )                                =                                                    Γ                (                                      ⁢                          I              ⁡                              (                u                )                                      ⁢                                              (        8        )            
Now the maximization is done on Q(I|I′)=E[ ln P(x|I)p(I)|y; I′]. Plugging the Gamma prior into Q, and solving for I(u), we get the following EM iteration for the maximum posterior estimation:
                                                                        I                ⁡                                  (                  u                  )                                            =                            ⁢                                                +                                                            ∑                      t                                        ⁢                                                                                  ⁢                                          E                      ⁡                                              [                                                                                                                                            x                                t                                                            ⁡                                                              (                                u                                )                                                                                      ❘                                                          y                              t                                                                                ;                                                      I                            ′                                                                          ]                                                                                                              +                                      ∑                                                                                  ⁢                                                                  ϕ                        t                                            ⁡                                              (                        u                        )                                                                                                                                                                    =                            ⁢                                                                    1                                                                  ⁢                                            +                                                                        ∑                          t                                                ⁢                                                                                                            ϕ                              t                                                        ⁡                                                          (                              u                              )                                                                                ⁢                                                                                                      [                                      +                                                                  ∑                        t                                            ⁢                                                                                          ⁢                                                                        y                          t                                                ⁢                                                                              ϕ                            ⁢                            ⁢                                                          (                              u                              )                                                        ⁢                                                                                          I                                ′                                                            ⁡                                                              (                                u                                )                                                                                                                                          ⁢                                                          ϕ                              t                                                        ⁢                            ⁢                                                          I                              ′                                                        ⁢                                                                                                                                ]                                ⁢                                  (                  10                  )                                                                                        (        9        )            
The EM update step can be formulated in matrix notation as follows. Let Φ be the matrix of the projections [φl(u)]t,u, and let I, I′, y, α and β be represented as column vectors. Eqn. 10 can be written in vector and matrix notations as:
                    I        =                              α            +                                          I                ′                            ·                              (                                                      Φ                    T                                    ⁢                                      y                                          Φ                      ⁢                                                                                          ⁢                                              I                        ′                                                                                            )                                                          β            +                                          Φ                T                            ⁢              1                                                          (        11        )            where the explicit multiplication and division denote element-wise operations, and where 1 is a vector (of the appropriate length) consisting solely of 1's.
Limited computational resources (i.e., when the entire projection matrix Φ cannot be kept in memory) may require breaking the update computation according to a partition of Φ into a set of sub-matrices (Φi). In that case the intensities can be updated gradually (using only one sub-matrix at each step) according to the following computation:
                    I        =                              α            +                                          I                ′                            ·                                                ∑                  i                                ⁢                                                                  ⁢                                                      Φ                    i                    T                                    ⁢                                                            y                      i                                                                                      Φ                        i                                            ⁢                                              I                        ′                                                                                                                                      β            +                                          ∑                i                            ⁢                                                          ⁢                                                Φ                  i                  T                                ⁢                1                                                                        (        12        )            where yi is the vector of observations that are obtained using the views of Φi.
In the context of image reconstruction from a constrained set of views, the utility of the reconstruction algorithms described above is limited. Many estimation algorithms, including EM and the Radon transform, require measurements from a complete set of views surrounding the imaged object. Although algorithms for EM with missing views have been developed, these algorithms are based on equally spaced views surrounding the imaged object, of which a number of views are not available. These algorithms do not provide a generalized solution for an unconstrained set of views, which may not be equally spaced or available for all directions surrounding the object.
Singular value decomposition (SVD) is a known technique for factorizing a rectangular real or complex matrix, with applications in signal processing and statistics. SVD may be considered a generalization of Eigenvalue decomposition to m*n matrices, whereas Eigenvalue decomposition is applicable only to square matrices.
SVD states that given the m-by-n matrix M whose entries are either from the field of real numbers or the field of complex numbers, there exists a factorization of the form:M=UDVT  (13)where U is an m-by-m unitary matrix, D is m-by-n with nonnegative numbers on the diagonal and zeros off the diagonal, and VT, the conjugate transpose of V, is an n-by-n unitary matrix. The elements along the diagonal of D are denoted the singular values. Such a factorization is called a singular-value decomposition of M. A common convention is to order the singular values Di,i in non-increasing fashion, so that the diagonal matrix D is uniquely determined by M.
The condition number of a matrix is defined as the ratio of the matrix's largest singular value to its smallest singular value. In numerical analysis, the condition number associated with a problem is a measure of that problem's amenability to digital computation. A problem with a low condition number is said to be well-conditioned, while a problem with a high condition number is said to be ill-conditioned. For example, the condition number associated with the linear equation x=My gives a bound on how accurate the solution y will be after approximate solution.
SVD may be employed for the solution of linear inverse problems. The inverse problem for a set of linear equations, x=My, is to calculate vector y from a known x and M. Using SVD, with M=UDV*, the problem may be restated as x=UDV*y. If D is an invertible matrix, y is obtained as:y=VD−1UTx  (14)where U and Vt are unitary matrices, and D is a diagonal matrix containing the singular values of M. Since U and Vt are easily transposable, a solution to Eqn. 13 is obtainable only if D is invertible. As D is a diagonal matrix of singular values e1 to en (where ek denotes the k-th singular value Dk,k), D−1 is a diagonal matrix of the reciprocals of the singular values. That is:
                              D                      -            1                          =                  [                                                                      e                  1                                      -                    1                                                                                                                                                                                                                                                                                                                                                                                                                                                e                  2                                      -                    1                                                                                                                                                                                                                                                                                                                                                                                                                              …                                                                                                                                                                                                                                                                                                                                                                                                              e                  n                                      -                    1                                                                                ]                                    (        15        )            
If the dimension of y is larger than the dimension of x, y cannot be determined directly but rather can only be estimated, which is often performed by an iterative procedure. If the condition number of D is very large, the multiplicative factors ei−1 will vary greatly, thus multiplying any measurement or calculation errors based on individual or linear combinations of elements of y, and decreasing the likelihood of convergence of the iterative estimation procedure.
Truncated SVD is a known technique for reducing sensitivity to inaccuracies or noise when solving a set of linear equations. In truncated SVD, the constraints associated with the smaller singular values are eliminated from the estimation. In terms of the inverse problem of Eqn. 14, this is accomplished for the i-th singular value by setting equal to zero in matrix D−1. Thus the lower valued singular values no longer affect the estimation process.
In order to achieve a reconstructed image which is adequate for medical diagnostic and treatment purposes, a reliable, high-resolution image of the tested object (i.e. body structure) must be obtained. Currently, reliable reconstruction algorithms are available only for complete data sets, which provide coverage of the entire volume. Such data is generally not available during medical imaging. Additionally, when high-resolution detecting units are used, their efficiency is relatively low, and the detecting units must remain at each position for a relatively long time in order to achieve a high probability of detection. Since during medical testing, measurements are generally performed at many locations as the detecting unit is moved relative to the observed body structure, the testing procedure generally requires a long time and is physically and emotionally difficult for the patient. Additionally, reconstruction is based upon a large quantity of data, and is a lengthy and computationally complex process.
There is thus a widely recognized need for, and it would be highly advantageous to have, an apparatus, system and method devoid of the above limitations.